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Abstract 

In this paper, we are concerned with improving the forecast capabilities of the Global approach 
to Time Series. We assume that the normal techniques of Global mapping are applied, the noise 
reduction is performed, etc. Then, using the mathematical foundations behind such approaches, 
we propose a method that, without a great computational cost, greatly increase the accuracy of 
the corresponding forecasting. 

I. INTRODUCTION Its relevance may be gauged by the existence 

of extensive studies in a great diversity of 
For any observed system, physical or oth- branches of knowledge, in physics as well as 
erwise, one generally wishes to make predic- in economics and the stock exchange, meteo- 
tions on its future evolution. Sometimes, rology, oceanography, medicine, etc. 
very little is known about the system. Pos- 

A time series is normally taken clS db set 

sibly, the dynamics behind the phenomenon 

of numbers that are the possible outcome of 

being studied is unknown, and one is given 

measurements of a given quantity, taken at 

just a time series of one (or a few) of its 

regular intervals. In reality, however, the as- 

parameters. Therefore, performing a time- 
sumption that the time series reflects in some 

series analysis is the best one can do in order 

way the underlying dynamics of the systems 

to learn the properties of the phenomenon. 

is worsened by the fact that the measured 
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may be due to a random external influence 
on a linear system, a noise (induced possibly 
by the measuring apparatus or other sources 
of contamination) which gets mixed with the 
desired information, thereby hiding it. But it 
may well be that they appear as a manifes- 
tation of low-dimensional deterministic chaos 
resulting from an intrinsic nonlinear dynam- 
ics governing the quantity under study (over 
which a random noise may also be superim- 
posed), with the characteristic sensitivity to 
initial conditions. 

If the time series is the only source of in- 
formation on the system, prediction of the fu- 
ture values of the series requires a modelling 
of the system's (perhaps nonlinear) dynami- 
cal law through a set of differential equations 
or through discrete maps. However, it is even 
possible that we do not know whether the 
measured quantity is the only relevant de- 
gree of freedom (frequently it is not) of the 
dynamical problem, nor how many of them 
there are. 

Both noise-contaminated linear and non- 
linear systems have nevertheless been stud- 
ied with success employing statistical tools, 
chaos-theory concepts, together with time- 
series analysis . Given a time series, one 
should ask first whether it represents a causal 
process or whether it is stochastic. Tools 
have been developed to decide upon this fun- 
damental question (the most common ones 

2 



are spectral analysis, Lyapunov characteris- 
tic exponents and correlation functions, see 
3], 0] ) . In the case of a series originated from 
a low- dimensionality chaotic dynamics, tra- 
ditional linear methods of analysis are not 
adequate, but an analysis apparatus was de- 
vised for applications to such nonlinear sys- 
tems I21 0| and we will not be concerned with 
stochastic processes in this paper. 

Methods for dealing with nonlinear time 
series fall mainly into two categories: local 
or global methods. Local methods are based 
on the assumption that, while in the long 
run nearby trajectories on the phase space di- 
verge considerably, they stay within the same 
neighborhood for a while. One may conjec- 
ture that to predict the next step in a time 
series, a good indication should come from 
the previous visits the system had made to 
the phase space neighborhood containing the 
"last point" of the series. An average of 
the behavior of the system for neighboring 
points, with a minimization of the distance 
in the phase space between them, gives good 
results for the next-step forecasting. 

Global methods, on the other hand, pos- 
tulate a functional form for the dynamics to 
be valid for any time. Usually one consid- 
ers polynomials of a suitable degree and one 
should devise a convenient way to estimate 
its coefficients. In this paper, we are going to 
concentrate in the global approach and, ac- 



tually, we will start from the global mapping properties of the corresponding attractor of 

itself, i.e., we are not going to be concerned the system we have to reconstruct it. 
with how the global mapping was generated 

(there are many standard approaches to do In l±A Takens used a method to recon- 

it) and we will not deal with noise reduc- struct the P hase s P ace - Vectors & ( with 

tion either (such considerations are impor- dimension "m") are reconstructed from the 



tant when determining the mappings, etc.). 
We will focus on a new method to, from any 
standard mapping one might have, improve 
the forecasting using it, without having to 
pay a very high computational price. 

Nonlinear analysis of Time Series relies 
not on the original maps of the dynamic 
system, but on its time- delay reconstruction. 
All discussions on the nonlinear treatment 
of Time Series make use of this reconstruc- 



Time Series Xi where Xi = x(U),i = 1,...,N 
as follows: 



6 = {x(ti),x(ti+p),...,x(ti+(m-l)p)} (1) 



tion scheme. There are already classical ref- 

Ojlsj]. This method allows one to reconstruct 
the phase space of the system with reasonable 
accuracy, using the information contained in 
the series only. 

Lorenz j^] showed that dynamic systems 
of low dimensionality could present strange 
attractors on their phase spaces. Takens [l0| 
proposed a method to reconstruct such phase 
spaces from the knowledge of a Time Series 



where m is the embedding dimension and p is 
the time lag (for definitions, see Based 
on the trajectories of the reconstructed at- 
tractor, we can study various topological in- 
variants of the system such as the Lyapunov 



ll|, etc. 



exponents, the generalized entropies 
We can also extract the underlying dynamics 
via a global modelling of the system. For 
example, one can try to obtain a low order 
Taylor series expansion for the system, thus 
obtaining a global mapping representing the 
system. We can use this mapping to perform 



a forecast of entries we ignore, i.e., in the fu- 

obtained from the system. He demonstrated i 

ture 1 . 

that the original attractor and the recon- 



structed one are characterized by the same 1 One can also do forecasting in a local version via 

analyzing the behavior of close vectors (to the one 



asymptotic 
acteristics 



properties and topological char- 



just before the one to be predicted) in order to 



11 1 . So, if we want to analyze the estimate the next (unknown) entry (see 0) 



II. AN ALGORITHM TO IMPROVE 
THE GLOBAL FORECASTING 

A. Stating the problem 

Suppose that the system can be modelled 
by a set of differential equations of low di- 
mensionality. What we would like to ob- 
tain is some kind of global map that, given 
any point of the phase space, could calculate 
a subsequent point of the trajectory. If we 
have known the set of differential equations 
(SDE) that models the system, we could find 
a solution (starting from an initial condition) 
by making a numerical integration through 
some map obtained from the SED (probably 
a Runge-Kutta map, a Taylor series one or 
an expansion in some function basis). For 
practical purposes (computers can not work 
with the infinity) a truncation must occur at 
some order of the series expansion. How- 
ever, if the truncation order is low, we can 
run away from the real solution in a few time 
steps (even if each time step is very small). 
For chaotic systems it is not used (in general) 
a Runge-Kutta expansion of degree less than 
four. This implies that the map generated 
present polynomials of high degree. Let's ex- 
emplify using one of the simplest chaotic sys- 
tem that exists, the Lorenz system: 

xi = cr(x 2 - £1), 



x 2 = —x 2 - Xi x 3 + Rxi, (2) 
x 3 = x 1 x 2 - bx 3 , 

where a, R and b are parameters and the sys- 
tem presents chaotic behavior for R > 24, 74. 

Why one of the simplest! Notice that 
this system possesses the minimum number 
of autonomous 2 differential equations perme- 
ating chaos: three 3 . Besides that, chaos is 
a phenomenon that only takes place in non- 
linear systems, and the smallest piece of non- 
linearity that we can add to a linear system 
in order to turn it non-linear is a quadratic 
term. 

Observe that the Lorenz system presents 
only two non-linear quadratic terms. Even 
in this simple case, as we will show, a Taylor 
series expansion of fourth order, will lead to 
a map of fifth degree in three variables. 

Consider the following initial condition: 

Xl (0) = x 10 , x 2 (0) = x 20 , x 3 (0) = x 30 . We 
can expand the corresponding solution as: 

(3) 

Since the system is defined by the equations 
^ = fi(x), we have that x * = fa = j { 4 , 

2 The time does not appear explicitly. 

3 In two dimensions we can not have chaos because 
the trajectories can not cross. 

A TTT1 • 

Where u represents — . 

at 



implying that: using the Lorenz system as a model for the 

Global Fitting scheme, let us suppose that 

d ( #A _dfi dft dx 3 _ ^ dfi 

dt \ ~d~t I ~ ~dt ~ £kc~~dt~ ~ ~dx we ve a *- ime Series produced from this 

(4) system (for instance, take one of the coordi- 

We can notice that, for the case of the nates of the system) After the usual phase 

Lorenz system, this process will increase the s P ace reconstruction ^ , say we want to have 

degree of the polynomials forming the map- a fourth order mapping (for the reconstructed 

ping by one for each order 5 . So, the map- system) with the same accuracy that could be 

ping corresponding to the forth order Tay- found on the fourth order Taylor expansion 

lor expansion is, at maximum, formed by for the Lorenz system. We would have to 

fifth degree terms. A polynomial mapping employ some minimization technique to de- 

of fifth degree implies a total of 168 coeffi- termine 168 coefficients. In practice, this is 

cients. Please remember that, as mentioned, a ve] T ni S n number making the whole proce- 

this is for one of the simplest chaotic dynamic dure computationally expensive. 

system cases (i.e., three-dimensional and only , r , , , , . c . 

v So, we are left with the hard choice of: ei- 

two non-linear (quadratic) terms). . , . , 

v ^ ' ' ther pay the computational price mentioned 

It is important to notice that, m Time Se- above and be very patient or try anc j decrease 
ries analysis, we do not have the dynamic the degree of the mapping . f course, there 
system to begin with. We, of course, will ig no such thing as a free meal The prke for 
consider that there is such a system behind the latter choice would be that the accuracy 
the series and we will look for determining would decrease (the corresponding Taylor ex- 
it. With the explanations above, we hope to pansion would be of lower order ). 
have made it clear that, even if the underly- 



Therefore, despite the fact that the global 
approach has many attractive features, such 



ing system is as simple as the Lorenz's one 

we will already have to face a great compu- 

. . i . i /• r , r , as the fact that, once it is determined it is ap 

tational task (it one wants to use forth order 



expansions - generally the minimum accuracy 
necessary for practical purposes) of determin- 
ing the 168 coefficients. With more detail, 



plicable to the whole series , one sees that the 
effective use of it can be difficult to achieve 
in practice. So, there is a clear demand for 
procedures that can, without increasing the 



Since the highest degree present in the functions 

/, g and h is quadratic, the derivatives (present on 6 In the case of Local mappings, we have to deter- 
I0J) are, at maximum, first degree polynomials. mine a mapping for each entry of the series. 
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degree of the global mapping, enhance the variables: 

accuracy of such mappings. In the next sub- x * = Fi(x,t), (5) 

section, before introducing one such attempt, 

where t is the group parameter. From Lie's 

we will talk about mappings. 

theory (13|, [lj, 



151 ] . we know that this group 



is the solution to a SDE defined by: 

B. Regarding Mappings 

Xi = fi(x), (6) 

In order to clarify the central idea of 
our proposed algorithm, let us make some where fi(x) = ^ \ t=0 and Xi = There- 
comments and present some results concern- fore, the transformation group (0) (i.e., the 
ing mappings representing the solutions for solution to the dynamic system ©) can be 
SDEs. obtained from the group generator defined as 

Consider the transformation group in n the operator X = YT%=i fi g|:> as follows: 



a.2 oo ±k 

x* = Fi(x, t)=x t + tX[ Xi ] + - } X 2 [x t } + ••• = £ TjX k [xi}. (7) 

- " k=0 



In this way, starting from a generic point In practice, the process of numerically solving 

Po, with corresponding coordinates X(p ), by the SDE can be summarized by choosing a 

choosing a time interval St, the transforma- small time interval (St <C 1) and truncating 

tion group © generates a mapping M that the series (jSJ) at some order n, thus obtaining 

takes a point on some given solution to the a mapping M given by: 
system and takes it to another such point 

that corresponds to a group parameter in- n ^k 

creased by St "•< P+ " = ^'^^ = JljT^M 

(9) 



x i(p+i) — Fi{ x (P),$t) — ^2 ~jT-^ k \- Xi (.P)\' where £j(p +1 ) approaches Xj(p +1 ) when St 



k=o fi/ - 

(8) 0. Defining the functions S k C{ clS 
6 



(et(£(p)) 

(6ei(x( P )) 
5 k ei(x {P) ) 



t 

5°£i(x(P))) = X i{P+1) -X i[P+1) = Tr ] Xk l X i(P) 

k=N+l K - 



= 5 1 e i (x( P ))^ = £j(f(p + i)) - £i(f(p)) 
= <J fc_1 e i (x( P+ i))-(J fc ~ 1 ei(f(p)), 



(10) 



where (k=2,. . . ), one can notice that 

= E *=< + o(^ 2 ) (ii) 

3=1 

and, generally, 



E 



5 fe+1 £i(f(P)) =5^(f ( p)) = 

6x { + 0(5xi 2 ). (12) 



dS*Ei(£(p)) 



dxi 



Since 5£ — > implies that 5xi —>■ 0, we can 
using (JT2J), enunciate the following result: 

5 k+1 e 



lim 

st^o S k e 



0. 



(13) 



where k is a positive integer. 

In the next subsection, based on this im- 
portant result, we will present an algorithm 
that enhances the predictive power of global 
mappings for Time Series. 

C. Mathematical Basis for the Algo- 
rithm 

Based on the above result (jT3j) , we have 
produced an algorithm that allows for im- 
proving the forecasting for the global fitting 
of a Time Series. 



As mentioned, we will suppose that the 
given Time Series is originated from phenom- 
ena that can be described by a low dimension 
dynamic system (Sp). After the phase space 
reconstruction [10j . we have a set of vectors 
defining a set of points along a single trajec- 
tory of the reconstructed systems (S r ) 7 . As 
usual, what we would like to determine is a 
global mapping M that would (with infinite 
precision) represent the solutions of the sys- 
tem S r . But, of course, in practice, what we 
can do is to produce a global mapping M 
through a procedure involving a minimiza- 
tion process 8 . If the Mapping M produces 
good forecasting for the series, that means 
that the coefficients present on M are close 
to the analogous ones present on the map- 
ping M which can be represented by the infi- 
nite series (JBJ) (and, ideally, it would describe 
S r with absolute precision). In that situa- 
tion, we would be in a similar position to the 

7 Takens ^3 has demonstrated that the system So 
and S r are tope-logically equivalent. 

8 In layman terms, what is done is to adjust the 
coefficients of the polynomial mapping (of a certain 
degree) to better reproduce the phase space points. 
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one presented on the last subsection (where 
we had just a truncated series because we 
knew the underlying SDE and could deter- 
mine the Taylor expansion). Why similar? 
In the "real" case we are dealing with now, 
we only have the series and have to deter- 
mine the mapping through a finite process 
and, therefore, the coefficients would not be 
exactly the same as in the truncated expan- 
sion of S r . So, defining functions and 5 k Si 
analogously to how we did in the last sub- 
section, we would expect that (JT3Jl would be 
valid. Actually, in the real world, the inequal- 
ity 

S k+1 e < 5 k e (14) 

is not valid for any positive integer k. The 
point is that, in actual calculations, 8t would 
be a finite value (not infinitesimal) At. So, 
at some integer value k, the inequality (pHj) 
would become 

A K+1 e « A K e. (15) 

The above reasoning allows us to build an 
easily applicable algorithm: Consider that 
we want to forecast the coordinate Xj (where 
% can take any value from 1 to the dimen- 
sionality of the reconstructed system) of a 
point P + 1 that immediately follows a given 
point P. In order to produce the mapping 
M, we use a certain number a + 1 of points 
that precede the point P + 1 (the points 
P, P — l, P — 2, . . . , P — a). Using this map- 
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ping, we can forecast the Xi coordinates for 
these a + 1 points. Let us call these a + 1 val- 
ues Wi. From these, we can define the func- 
tions A k Ei (analogously to the functions (fTUJl 
in subsection IIIBjl . 

A°£i(x {J) ) = x i{J) -x i{J) 

A 1 e i (x( J )) = A°£i(f ( j)) - A°e i (f( J _i)) 

A k E l (x {J) ) EE A fe - 1 £ i (f (J) )-A fe - 1 £ i (x (J _ 1) ), 

: i (16) 

where (k = 0, . . .) and (J = P — a + k, . . . , P). 
Using these definitions, we can determine the 
values for k where we have A k+1 e ~ A k e 9 
and, using this knowledge, we will see that 
we can improve the forecasting generated by 
the mapping M. Let us clarify what we 
mean: if we want to forecast the value for 
the coordinate X{ of the point P + 1, we 
may use the global mapping M that would 
produce the forecast Xjrp+i). We know that 
x i(P+i) - %i(P+i) = A°£j(f( P+1 )) and, there- 
fore, 

%i(p+i) = z*(p+i) + A°e i (f( P+ i)). (17) 

Notice that we do not know the value 
for A°£j(5f(p + i)). But we know that 
A 1 e < (f(p +1 )) = A £ i (f(p + i))-A°£ i (f(p)), im- 

9 There is a finite range for the values for k in which 
that happens. After a certain value, the A fe ~ 1 £i 
start to diverge. 



plying that 

A°Ei(x( P+1 )) = A Ei(x( P )) + A l Ei(x (P+1 )). 

(18) 

Let us examine this: we know the value 
for A £j(x(p)) (i.e., £j(p) — £j(p)) but we 
do not know A 1 £j(x(p + i)). However, if (P) 
and (P + 1) are sufficiently close (such that 
A 1 ^ <C A°£j), we can expect that we will 
gain information when substituting (jl8j) into 
(fT7j) obtaining 

x i{P+1) = x i{ p +1) +A°£: i (f(p)) + A 1 £: i (x(p + i)). 

(19) 

Why do we gain information? If we 
compare (fTTj) to (fT9*j) . we can observe that 



x i{P+ i) = x i( p + i) + A°£j(f ( p)) + A^j^p) 

Therefore, we can build a simple algorithm 
to improve the prediction Xj(p +1 ), obtained 
with mapping M: we determine the integer 
k for which the approximation starts to fail, 
i.e., A k+1 Ei ~ A h Ei, then we neglect the term 
A fc+1 £ i (a ? (p + i)) and end up with 

#i(p+l) = Xj( P+1 ) + A £j(x(p)) + 

A 1 e i (i> ) ) + -- + A fc £ i (x> ) ). (22) 

The remaining question is: How to de- 
fine A k+1 £i ps A k e{! Let us elaborate the 
analysis just made above. We are interested 



the unknown term in (fTTj) is A°£j(x(p + i)) 
which is (by hypothesis) much bigger than 
the unknown term in (fTT?|) : A 1 £ i (x(p + i)). 
So, the term A°£j(a;(p)) is a correction to 
aJj(P +1 ). Analogously, we have A 2 £j(x(p + i)) = 
A^j^p+i)) - A 1 £ i (x( P )), implying that: 

A 1 £ i ( y X (P+1) ) = A 1 £ i (x {P) ) + A 2 £i(f ( p +1) ), 

(20) 

substituting this into (fT9*|) . if (P) and (P + 
1) are sufficiently close such that A 2 Ei 
A 1 ^^, we would have a second order correc- 
tion to x~ifp + xy Actually, when the relation 
A k+1 Ei <^ A k Ei applies, we can further cor- 
rect £i( P+ i), i.e., 

+ • • • + A k E, l {x {P) ) + A fe+1 £,(f ( p +1) ). (21) 
I 

in using an approximation, a kind of Tay- 
lor series expansion, when trying to forecast 
the Time Series, what one might expect from 
such a situation? In a perfect world, the 
terms in the series would, gradually, become 
smaller in an infinite fashion. Of course, as 
already mentioned above, we are dealing with 
a real series, where each entry is not infinites- 
imally apart the previous one and is, actually, 
finitely separated. How "finitely separated" 
depends on the particular series under study 
and, being more rigorous, on the particular 



section of the series we are considering. This 
translates to the fact that, if one considers the 
absolute values of the differences A k Ei, they 
will decrease with increasing values for k un- 
til this value reaches the magnitude defined 
by the "non- infinitesimal" character of the 
Time Series we have just emphasized, where 
this character will then make the values for 
the differences oscillate (for a while) around 
this magnitude (since this magnitude would 
dominate over the initial tendency of the dif- 
ferences to decrease). With the increasing 
values for k, this initial tendency of the dif- 
ferences to decrease will cease as our approx- 
imation (Taylor like) stars to diverge from 
the actual value for the series. We will then 
see the absolute values for the following dif- 
ferences start to increase and rapidly diverge. 
That clearly, if one thinks in plotting the (ab- 
solute) values for the differences, defines a 
plateau where A h+1 6i ~ A k 6i and our above 
introduced method will work at its best. 

D. The steps of the algorithm 

Consider that we have already recon- 
structed the phase space from the Time Se- 
ries under study and that we want to forecast 
the P + 1 entry (P is the last known value of 
the series). This entry corresponds to a coor- 
dinate of a reconstructed vector on the phase 
space (as usual). Using a global mapping M, 



obtained via standard k-fold validation pro- 

n 

cedures |16|, we do the following: 

1. Set n = 10. 

2. We calculate the absolute value for the 
functions A k Ei (see eq. (JTHJ)) up to k = n 
for the point P. 

3. We check to see if we have already 
found the plateau, i.e., we look for the 
value of k for which |A fc £j| < |A fc+1 £j|. 
Please note that this checking can be 
very easily automatized. 

4. If the checking returns false we set 
n=n+10 and return to step 2. Other- 
wise we would have found the corrected 
value for Xj(p +1 ) as: 

x i{P+l) = X i(P+1) + A°£ i (x (P) ) + 

A 1 s l (x iP) ) + --- + A k e i (x {P) ). (23) 

III. APPLICATIONS 

In this section, we are going to present 
two applications of the above introduced 
improved forecast method. We will start 
by introducing the Time Series in question, 
present the reconstruction parameters and 
the associated Global Mapping. We then will 
proceed to the algorithm, following the steps 
just introduced and compare the average per- 
formance for the "usual" and the improved 
approaches. 
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A. Application 1: Lorenz 

1. The Time Series 

This is an "academic" application in the 
sense that it is, certainly, originated from a 
dynamic system and we actually even know 
which one. But it is important in order for 
us to see the ideas of the improved method 
working on an arena that suits it very nicely. 

The Time Series was generated taking the 
consecutive values for the X\ coordinate of the 
Lorenz system (see eq. ©), starting from the 
initial condition x\ = —0.3336666667, X2 = 
-0.3336666667, x 3o = 21.9996666667, us- 
ing an eighth order Runge-Kutta numerical 
integration The Series presents 600 en- 
tries (please see figure ((TJ) for a plotting of 
this Time Series). 

Now, in order to apply the Global Anal- 
ysis (jll 0) to this Time Series, we have to 
reconstruct the phase space. To do that, we 
need to determine the relevant parameters, 
namely the time-lag and the embedding di- 
mension (please see [3]). For this present 
case, the reconstruction parameters are time- 
lag = 6 and embedding dimension = 3. So, in 
the remaining of this subsection, we will call 
these three dimensions of the reconstructed 
phase space for the Lorenz system (x,y,z). 
In real life, we use the whole Time Series we 
know/measure to produce the Global Map- 

11 



ping and use it to predict future (unknown) 
entries. Here, in order to evaluate the ac- 
curacy of the predictions we obtain using a 
regular Global Fitting and our Improved one, 
we are going to use an initial portion of the 
Series to generate the Mapping and the other 
(remaining) portion of the Series as our test- 
ing ground, i.e., we will apply our mappings 
to entries in that region and compare it to 
the actual values to see how the mappings 
fared. In the present case, the first 140 en- 
tries constitute our portion of the Series used 
to build the Mapping up. Basically, we use 
all the vector reconstructed from these en- 
tries and produce a quadratic fitting minimiz- 
ing the distances from this fitting (when ap- 
plied to each vector) to the actual values via, 
for instance, a least mean square procedure. 
Actually, we also have used an improvement 
(a very standard one) called a k-validation. 
In layman's language, basically what this k- 
validation does is to average up several map- 
pings. Doing all this, the global mapping we 
have derived (and to be used on this applica- 
tion henceforth) is: 

M = 1.317833301 x - 0.005266089766 a; 2 
+0.07580676400 xy - 0.1245478927x2 - 
0.01839588238 y 2 + 0.06578287850 yz - 
0.01025562766 z 2 - 0.4700554502 y + 

0.1415056465 z. 

(24) 



2. The inner works of the improved forecast 
algorithm 

Let us now, using two generic points from 
the Series, exemplify the workings of our im- 
proved method. 

Consider the entries P = 316 and P = 
533, with respective values of —1.370578116 
and 6.860383245. The values for the entries 
P = 317 and P = 534, the "next" entry for 
each case considered here, are —1.041455029 
and 7.225654731. Let us see how the "usual" 
Global Fitting fares in these entries. Us- 
ing the mapping presented on (f53j). we get 
the following forecasting for the entries P = 
317 and P = 534: -0.782644049 and 
7.062374264. These present a "percentage 
error" (given by \{value — forecast) /value\) 
of 24.85090309 and 2.259732482 respectively. 
How about the improved method? 

In order to apply our method we have to 
find the plateau by finding the value for k to 
which |A fc £:j| < |A fc+1 £j|. Let us do that for 
the couple of points chosen above: 

• P=316 

As "prescribed" above, what we have 
to do is, by looking at table (JTJ), second 
column, determine at which value of k 
A k e(x(p-)) stops decreasing for the first 
time (and begin the oscillations we have 
mentioned in section Hi C|) . From table 
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(0), we see that happens for k = 5. Us- 
ing this into equation (J2HJ), we find (see 
table (jTJ)) that the "percentage error" 
for our method is 0.0004798095. 

• P=533 Again, what we have to do is, 
by looking at table (fTTJ) , second col- 
umn, determine to which value of k 
A fc e(x(p)) stops decreasing for the first 
time (and begin the oscillations we have 
mentioned in section III C|) . From table 
(jTTJ), we see that happens for k = 3. Us- 
ing this into equation (J23|) , we find (see 
table (jTTJ) ) that the "percentage error" 
for our method is 0.0001986533. 

As we have mentioned in section pi C|) . 
we expect the absolute values of A k e(x(p)) to 
oscillate when \A k Ei\ ss \A k+1 Ei\. That fact 
is illustrated, for the entries P = 316 and 
P = 533 respectively, on figures (j2J) and (jHJ). 

3. Performance Comparison 

The reader may ask: why these two en- 
tries above? Fair enough, they are not special 
at all. So, in order to confirm the fact that 
our new approach may be an advantage, let 
us make a general survey of the entries on 
the Time Series. We take 21 entries, equally 
distributed, on the last part (not used when 
producing the Global Mapping) of the Time 



Series. The results are presented on table 

The idea behind of presenting the results 
for points equally spaced on the entire Time 
Series (meaning the entire testing ground 
defined above) was to provide the informa- 
tion on all the Time Series, i.e., it is very 
important (for many Series) the section in 
which the analysis is carried out. So, we 
have decided to present the results for many 
points, evenly distributed along every sec- 
tion of the Time Series. But, for complete- 
ness, we will present the average percentage 
error (for the improved Global fitting) for 
the whole testing ground for the Time Se- 
ries and for the 21 entries used on table 
The percentage error for the whole series is 
.1601961683e — 1 and for the 21 entries on 
the table is .3385849199e - 1. Both are com- 
patible, showing that the chosen 21 are rep- 
resentative of the totality of the possibilities. 
The percentage error for the "regular" global 
fitting is (for the whole series) 10.58939930. 

As can be seen, our method is a great im- 
provement of accuracy when compared with 
the "plain" Global Fitting. To help in this 
analysis we present figure (j3J) where we plot 
\n(A GF /A IGF ), where A GF and A IGF are, re- 
spectively the percentage errors in the Global 
Fitting and the Improved Global Fitting. As 
can be seen from the figure, most of the IGF 
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errors are smaller than e 4 times the GF er- 
rors. 



B. Application 2: Heart beat 

1. The Time Series 

Let us now deal with a more "real" exam- 
ple, where we deal with data extracted from 
Nature, we do not know the system behind 
the phenomenon, etc. The following Time 
Series was obtained 10 from measurements of 
the heart beat rate in a person performing 
many different activities. The Series presents 
1744 entries (please see figure (0) for a plot- 
ting of this Time Series). 

In order to produce the Global Mapping 
for this case, we have proceeded in the same 
fashion as we did in the Lorenz System Time 
Series application. So, we will not repeat the 
whole explanation of the procedures involved 
here. Please refer to section IV-A above. For 
this application, the reconstruction parame- 
ters are time-lag = 10 and embedding dimen- 
sion = 3. So, in the remaining of this subsec- 
tion, we will call these three dimensions of the 
reconstructed phase space for the Heart beat 
data (x,y,z). The global mapping we have 
derived (and to be used on this application 

10 http://ecg.mit.edu/time-series/ 



henceforth) is: 



the couple of points chosen above: 



M= -1.172534275 z- 0.2617292220 z 2 + 
0.4661889468 yz - 0.3426537822 y 2 - 
0.1107944588 xz + 0.3574412776 xy - 
0.1136908621 x 2 + 15.65933419 z- 
12.98866231?/. 

(25) 

2. The inner works of the improved forecast 
algorithm 

Let us now, using two generic points from 
the Series, exemplify the workings of our im- 
proved method. 

As in the previous application, consider 
the entries P = 737 and P = 1016, with re- 
spective values of 89.18875624 e 94.25098125. 
The values for the entries P = 738 and P = 
1017, the "next" entry for each case consid- 
ered here, are 89.16743126 and 94.28981563. 
Let us see how the "usual" Global Fit- 
ting fares in these entries. Using the map- 
ping presented on (J25J) . we get the follow- 
ing forecasting for the entries P = 738 and 
P = 1017: 90.996977 and 82.611004. These 
present a "percentage error" (defined above) 
of 2.051809404 and 12.38607961 respectively. 
How about the improved method? 

In order to apply our method we have to 
find the plateau by finding the value for k to 
which |A fc £j| < |A fc+1 £j|. Let us do that for 
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• P=737 

As "prescribed" above, what we have to 
do is, by looking at table (jIV|) . second 
column, determine to which value of k 
A fe e(x(p)) stops decreasing for the first 
time (and begin the oscillations we have 
mentioned in section III Cj) . From table 
PV|) . we see that happens for k — 1. 
Using this into equation (J2HJ), we find 
(see table (|IV)) ) that the "percentage 
error" for our method is 0.3970473916. 

• P=1016 

Again, what we have to do is, by look- 
ing at table (jVj), second column, de- 
termine to which value of k A k e{x^p)) 
stops decreasing for the first time (and 
begin the oscillations we have men- 
tioned in section ITT"C]) . From table (JVJ), 
we see that happens for k = 3. Using 
this into equation (j2Hj) . we find (see ta- 
ble (|Vj)) that the "percentage error" for 
our method is 1.524959626. 

As in the previous application, we expect 
the absolute values of A k s(x(p)) to oscillate 
when |A fc £j| « \A k+1 Ei\. That fact is illus- 
trated, for the entries P = 737 and P = 1016 
respectively, on figures © and ((7|). 



3. Performance Comparison 



IV. CONCLUSION 



Let us make the general survey of the en- 
tries on this Time Series. We take 26 entries, 
equally distributed, on the last part (not used 
when producing the Global Mapping) of the 
Time Series. The results are presented on 
table (jvTj) . 

The idea behind of presenting the results 
for points equally spaced on the entire Time 
Series (meaning the entire testing ground 
defined above) is the same one explained 
on the section regarding the Lorenz System. 
The percentage error for the whole series is 
1.877994467 and for the 26 entries on the 
table is 1.429515613. Both are compatible, 
showing that the chosen 26 are representa- 
tive of the totality of the possibilities. The 
percentage error for the "regular" global fit- 
ting (for the whole series) is 5.971546764. 

As can be seen, in the majority of cases, 
our method is, for this more "realistic" case, 
also a great improvement of accuracy when 
compared with the "plain" Global Fitting. 
To help in this analysis we present figure (jHJ) 
where we plot ^(Ag^/A/gf), where Aqf 
and AfQf are, respectively the percentage er- 
rors in the Global Fitting and the Improved 
Global Fitting. As can be seen from the fig- 
ure, most of the IGF errors are more than 
five times smaller than the GF errors. 



There is a huge demand for improving 
methods that do not cost too high a com- 
putational price to achieve desired levels of 
accuracy in Time Series Analysis. 

Here, we have presented one such method. 
The basic rational behind it is that we can 
make use, as explained in section |H1 of the 
underlying (assumed) low- dimensionality dy- 
namics to correct our forecast. It is impor- 
tant to mention that, in order to apply the 
method, one does not have to quantify the 
hyperbolicity (or the low- dimensionality, for 
that matter) of the Time Series. The steps of 
the procedure will take (automatically) care 
of stopping when this hyperbolicity "spoils" 
the correcting power of the method. So, 
the algorithm is secure. It is also useful to 
remember that our efforts here are aimed 
to avoid the computational cost of the fit- 
ting/minimizing procedures. So, our method 
is not equivalent to fittings, with the same 
computational cost, in any shape or form. 

We have presented two applications of our 
method: The first one is a (we are going 
to call it) pure low dimensional known sys- 
tem, from where we generated a Time Series. 
The reason for this application is to use the 
method on a controlled arena, i.e., we can 
see the method working at its best. What 
do we mean by its best? Could not have we 
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gotten better results than the ones presented 
in section MI AP Of course we could have, 
for instance, if we have made the Time Series 
more "dense", i.e. if we have used smaller 
values for At, of course, the results would 
be better. Indeed, we can do the same in- 
definitely up to infinite precision. What we 
mean by "its best" is the fact that there is 
not, for sure, any high dimensional behavior. 
We have then demonstrated that the ideas 
behind our method work quite nicely. 

The second application corresponds to a 
Time Series obtained from measurements, 
i.e., we do not have any prior knowledge 
about the (possible) dynamic system under- 
lying it. We have found that, after the usual 
techniques have been used to produce the 
Global mapping, we could improve the fore- 
cast capabilities of the fitting quite a bit (see 
section Ull Bj) . thus demonstrating the practi- 
cality of our approach on a uncontrolled sit- 
uation. 

Our method has, of course, its limitations. 



Perhaps the most obvious one is the fact that 
it won't help much in the case where the Time 
Series is "sparse", i.e., as we have mentioned 
just above, as At becomes large, the method 
won't work. The limitation so far is that we 
do not have a criteria, as yet, to, just by 
quickly inspecting the Time Series, determine 
if our method applies well or not. One has to 
have a go and, in a testing arena, verify if the 
method is improving things. 

That leads to future work: produce a fast 
algorithm to test the time series for appli- 
cability (or not) of the method. One other 
possible line of research to be pursued is to 
improve our algorithm in the sense of using 
more information contained on the plateau 
than we are using now. So far, we are tak- 
ing the first piece of data on the plateau but, 
as we have explained in section Hi CI the val- 
ues for the corrections will oscillate from that 
point on. It is reasonable to look for an algo- 
rithm to extract information from this oscil- 
lation. 
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k 


\A k e(x (P) )\ 


IGF 


error 


1 


0.015127123 


- 1.042168004 


0.06845950907 


2 


0.000902659 


- 1.041265345 


0.01821336445 


3 


0.000202785 


- 1.041468130 


0.001257951581 


4 


0.000032905 


- 1.041435225 


0.001901570346 


5 


0.000014807 


- 1.041450032 


0.0004798094839 


6 


0.000018502 


- 1.041468534 


0.001296743462 


7 


0.000011662 


- 1.041456872 


0.0001769639541 


8 


0.000009405 


- 1.041447467 


0.0007260995232 


9 


0.000006933 


- 1.041454400 


0.00006039627084 


10 


0.000006450 


- 1.041460850 


0.0005589295589 


11 


0.000005072 


- 1.041455778 


0.00007191861186 


12 


0.000011998 


- 1.041443780 


0.001080123451 


13 


0.000020272 


- 1.041423508 


0.003026630927 


14 


0.000055212 


- 1.041368296 


0.008328060030 


15 


0.000149033 


- 1.041219263 


0.02263813544 


16 


0.000346170 


- 1.040873093 


0.05587720869 


17 


0.000723317 


- 1.040149776 


0.1253297515 


18 


0.001398781 


- 1.038750995 


0.2596400156 


19 


0.002499815 


- 1.036251180 


0.4996710232 


20 


0.004042167 


- 1.032209013 


0.8877979118 



TABLE I: In this table, we plot the \A k e\, for 
the entry 316, for the Lorenz System Time Se- 
ries. IGF is our improved global fitting result 
corresponding to the particular value of k. 
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k 


A k e(x (P) ) 


IGF 


error 


1 


0.005833355 


7.225283437 


0.005138551644 


2 


0.000348756 


7.225632193 


0.0003119163708 


3 


0.000008184 


7.225640377 


0.0001986532783 


4 


0.000008912 


7.225649289 


0.00007531497425 


5 


0.000003874 


7.225653163 


0.00002170045565 


6 


0.000001257 


7.225654420 


0.000004304108231 


7 


0.000000331 


7.225654751 


0.0000002767915261 


8 


0.000000150 


7.225654901 


0.000002352727972 


9 


0.000000251 


7.225655152 


0.000005826461624 


10 


0.000000532 


7.225655684 


0.00001318911622 


11 


0.000001007 


7.225656691 


0.00002712556956 


12 


0.000001712 


7.225658403 


0.00005081892419 


13 


0.000002692 


7.225661095 


0.00008807506360 


14 


0.000004115 


7.225665210 


0.0001450249201 


15 


0.000006683 


7.225671893 


0.0002375148085 


16 


0.000012706 


7.225684599 


0.0004133604651 


17 


0.000028487 


7.225713086 


0.0008076084753 


18 


0.000068998 


7.225782084 


0.001762511561 


19 


0.000166009 


7.225948093 


0.004060005784 


20 


0.000380304 


7.226328397 


0.009323252011 



TABLE II: In this table, we plot the |A fc e|, for 
the entry 533, for the Lorenz System Time Se- 
ries. IGF is our improved global fitting result 
corresponding to the particular value of k. 
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N 


GF error 


IGF error 


300 


2.949988321 


0.0007879476836 


310 


9.384955078 


0.003895726200 


320 


602.9300615 


0.6329500170 


330 


4.588145520 


0.00003360765471 


340 


0.6804060144 


0.000006172635711 


350 


1.799833141 


0.000002478938512 


360 


1.908426262 


0.00003660198120 


370 


1.992750955 


0.00003095127924 


380 


5.351353323 


0.005676508421 


390 


15.14897504 


0.00006055593702 


400 


17.40670926 


0.03531034693 


410 


11.79516640 


0.00003410572689 


420 


6.417051601 


0.000009213921336 


430 


3.561881518 


0.0000003239205212 


440 


2.472699706 


0.0000008353740914 


450 


2.070219097 


0.000002536966462 


460 


2.807205469 


0.00004274128369 


470 


9.466251097 


0.03219442293 


480 


17.78540068 


0.00003404766212 


490 


15.54182977 


0.00001293826149 


500 


9.424552928 


0.0000008546217942 



TABLE III: Comparison between the Global Fit- 
ting and the Improved Global Fitting for the 
Lorenz Time Series 
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k 


A k e(x (P) ) 


IGF 


error 





1.40709476 


89.58988224 


0.4737727375 


1 


0.06841402 


89.52146822 


0.3970473916 


2 


0.11034922 


89.63181744 


0.5208024650 


3 


0.44228546 


90.07410290 


1.016819288 


4 


0.56575633 


90.63985923 


1.651306928 


5 


0.43932335 


91.07918258 


2.144001787 


6 


0.00174856 


91.07743402 


2.142040802 


7 


0.99434282 


90.08309120 


1.026899538 


8 


3.10231730 


86.98077390 


2.452304983 


9 


7.64654504 


79.33422886 


11.02779598 


10 


17.67633668 


61.65789218 


30.85155498 



TABLE IV: In this table, we plot the |A fc e|, for 
the entry 737, for the Time Series with the Heart 
Beat data 




200 400 600 



FIG. 1: Lorenz Time series. The horizontal axis 
marks the position of the entry (i) and the ver- 
tical on the value for the entry (X(i)) 
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k 


A k e(x (P) ) 


IGF 


error 





9.89683125 


92.50783525 


1.889896982 


1 


2.18086125 


94.68869650 


0.4230370664 


2 


1.05010575 


95.73880225 


1.536737144 


3 


0.01110500 


95.72769725 


1.524959626 


4 


0.49314562 


95.23455163 


1.001949143 


5 


0.37340972 


94.86114191 


0.6059257579 


6 


0.43065248 


95.29179439 


1.062658521 


7 


2.97396166 


98.26575605 


4.216723082 


8 


10.77153154 


109.0372876 


15.64057780 


9 


32.44747527 


141.4847629 


50.05306984 


10 


86.66665506 


228.1514179 


141.9682512 



TABLE V: In this table, we plot the |A fc e|, for 
the entry 1016, for the Time Series with the 
Heart Beat data 
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N 


GF error 


IGF error 


500 


3.397125034 


2.168865202 


540 


9.089213645 


1.889289807 


580 


5.874155830 


0.1858455792 


660 


1.844997944 


1.332964363 


700 


0.3583492837 


0.01355789862 


740 


2.400959523 


0.1672692403 


780 


0.4937771060 


0.1540773590 


820 


1.249041343 


0.1600760146 


860 


10.26039663 


0.8787867282 


900 


6.962087707 


0.8589041386 


980 


4.565823761 


0.2342961164 


1020 


14.67475919 


3.116867623 


1060 


1.249489572 


0.5393005133 


1100 


1.075608307 


2.987879519 


1140 


0.2169661915 


0.08177074130 


1180 


13.66135448 


0.4848516873 


1220 


2.732512685 


0.8664583752 


1260 


6.861097043 


0.6660946812 


1340 


9.400776847 


0.9535935739 


1380 


1.392042607 


0.2352261816 


1420 


0.8176069275 


0.5917808721 


1460 


2.813583072 


0.2695624030 


1500 


6.384669758 


6.651398095 



TABLE VI: Comparison between the Global Fit- 
ting and the Improved Global Fitting for the 
Time Series with Heart Beat data 
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n 1 1 1 1 1 r~ 

3 K 7 9 11 13 15 
K 



FIG. 2: The plot shows the values of the \A k e\ 
against the number k, for the entry 316, for 
the Lorenz System Time Series. In the x-axis, 
marked with the letter K, is the value of k that 
our procedure defines as the beginning of the 
plateau. 
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K 

FIG. 3: The plot shows the values of the \A k e\ 
against the value of k, for the entry 533, for 
the Lorenz System Time Series. In the x-axis, 
marked with the letter K, is the value of k that 
our procedure defines as the beginning of the 
plateau. 
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FIG. 4: The plot is for \n(A GF /A IGF ) against 
the position in the Time Series (i). The line 
marks the threshold where, above it, Ajqf 
starts to be smaller than e -4 times Aqf- 
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FIG. 5: Heartbeat data. The horizontal axis 
marks the position of the entry (i) and the ver- 
tical on the value for the entry (X(i)) 
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FIG. 6: The plot shows the values of the \A k e\ 
against the value of k, for the entry 737, for the 
Heartbeat Time Series. In the x-axis, marked 
with the letter K, is the value of k that our pro- 
cedure defines as the beginning of the plateau. 
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FIG. 7: The plot shows the values of the \A k e\ 
against the vale for k, for the entry 1016, for the 
Heartbeat Time Series. In the x-axis, marked 
with the letter K, is the value of k that our pro- 
cedure defines as the beginning of the plateau. 
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FIG. 8: The plot is for ln( Aqf / ^igf) against 
the position in the Time Series (i). The solid 
line marks the threshold where, above it, Ajgf 
starts to be the half of Agf and the dotted line 
the threshold where, above it, Ajgf starts to be 
be the fifth of Aqf. 



28 



